Microrheology of supercooled liquids in terms of a continuous time random walk 
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Molecular dynamics simulations of a glass-forming model system are performed under 
application of a microrheological perturbation on a tagged particle. The trajectory of 
that particle is studied in its underlying potential energy landscape. Discretization of 
the configuration space is achieved via a metabasin analysis. The linear and nonlinear 
responses of drift and diffusive behavior can be interpreted and analyzed in terms 
of a continuous time random walk. In this way the physical origin of linear and 
nonlinear response can be identified. Critical forces are determined and compared 
with predictions from literature. 
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I. INTRODUCTION 



In microrheological experiments a physical system is brought out of equilibrium by a 
rather simple perturbation which is acting on a small fraction of its constituents^. These 
experiments are of specific interest because the application of a well-defined perturbation 
allows, first, to give predictions about the non-equilibrium behavior of the system and, 
second, to gain additional information about the equilibrium properties of the system by 
analyzing the transition to the non-equilibrium state. 

In supercooled liquids in its (quasi-)equilibrium state one can observe a lot of unique 
dynamical properties like non-exponential relaxation or violation of the Stokes-Einstein re- 
lation^. Typically this is interpreted in terms of the presence of dynamic heterogeneities; see^ 
and references therein. In recent years also the microrheological properties of supercooled 
liquids have been studied experimentally^ as well as theoretically The occurrence of a 
nonlinear dynamical response was reported for different dynamical quantities. As expected 
for general grounds^, the nonlinear response strongly increases at lower temperatures. 

For the equilibrium dynamics of supercooled liquids it has been shown that a descrip- 
tion of the dynamics in terms of transitions between metabasins (MB) yields important 
additional information about the nature of the slow dynamics at low temperatures^^. It 
has been shown, that this discrete dynamics fulfills all criteria of a continuous time ran- 
dom walk (CTRW ) 13 * 14 ! if applied to the description of a small system (O(10 2 ) particles). 
Most importantly, it turns out that the diffusivity of a small system only shows very small 
finite-size effects^ whereas the structural relaxation time displays large effects. This can be 
interpreted via facilitation effects^. Going down in temperature may in principle modify 
the spatial and the temporal aspects of the CTRW. Interestingly, it turns out that the slow- 
ing down at low temperatures is exclusively determined by the temporal contributions^. 
Furthermore one can, e.g., express the complete wave vector dependence of the structural 
relaxation part of the incoherent scattering function in terms of the properties of the waiting 
time distribution, characterizing the CTRW^. Thus, the mapping of the MB dynamics on 
the CTRW is a powerful concept for the description of supercooled liquids. Unfortunately, 
the latter approach can be only used for small systems^. The qualitative reason is that 
very large systems can be decomposed into basically independent smaller systems^. Thus, 
spatial information about the dynamic processes becomes essential which, however, is not 
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available in configuration space. More formal arguments can be found irP. 

On a qualitative level the external force gives rise to a tilting of the potential energy 
landscape for the tagged particle. The key goal of the present work is to show how the 
CTRW properties are modified upon this tilting. This allows us to characterize the onset of 
non-linearity in the mobility of the tracer particles in terms of specific CTRW properties. 
Furthermore, it will be demonstrated that a small system, analogous to the equilibrium 
case, only displays small finite size effects, even in a high nonlinear regime. Therefore, a 
closer understanding of microrheological effects in small systems allows us to unravel the 
underlying physics of large systems as well. 

The paper is organized as follows: In Sect. [IT] we describe details of the simulation. 
Sect. Hill contains the results and their discussion. We conclude in Sect. [TV] 

II. SIMULATIONS 

We present the results of molecular dynamics (MD) simulations of a binary Lennard- 
Jones mixture, which consists of two types of particles, A and B, with an A:B ratio of 
80 : 20. This system is known to be a prototype of a glass-forming system 16 . To be able 
to simulate system sizes as small as 65 particles the cutoff radius is reduced to 1.8 rather 
than 2.5 (in dimensionless LJ-units) 17 . A microrheological perturbation is introduced to the 
system by randomly selecting one A type particle and pulling it with a constant force F 
along a certain direction. To ensure that the system resides in a stationary state, the system 
is equilibrated under application of the force to the tracer particle. A constant temperature 
during our simulations is achieved by coupling the system to a Nose-Hoover thermostat?^. 

A key element of our analysis is the tracking of minima of the PEL which the system 
explores during its time evolution and which displays the impact of the external force. 
We access these minima by minimizing the particles coordinates in certain time intervals 
with respect to its potential energy Furthermore, following the procedure as outlined 
in RefP^ the minima are grouped together to metabasins (MB). As mentioned above the 
potential energy is minimized for the complete potential energy, containing also the effect 
of the external force. In practice it turns out that for somewhat larger forces (F > 10) the 
minimization routine becomes instable. In any event, no larger forces are needed since the 
onset of non-linear effects occurs at much smaller values of F . 
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In MD simulations one generally prefers to simulate systems as large as possible to avoid 
any unphysical behavior due to the finite-size of the system. However, since the CTRW 
approach in configuration space is only valid for small system sizes^ we focus on a small 
system with N = 65 (BLJM65) which is close to the smallest system size without significant 
finite-size effects for the diffusivity 19 . Due to linear response theory this automatically 
implies that also the linear regime of the mobility in the microrheological setup should not 
display relevant finite-size effects. Whether or not this holds also for the non-linear regime 
will be checked by comparison with the simulations of a larger system with 1560 particles 
and a geometry of 3 : 1 : 1, the long side pointing along the force direction (BLJM1560). 

III. RESULTS AND DISCUSSION 

A. Drift velocity in real space 

One central quantity of the analysis of microrheological properties is the drift velocity v 
the tracer particle obtains due to the interplay of permanent acceleration along the force 
direction and friction effects. Formally, the drift velocity is given by the stationary long-time 
limit of the displacement of the tracer particle parallel to the force direction x\\(t) relative 
to the time. 

v = hm X " (t) - X " (t0) . (1) 

t->oo t — to 

For BLJM65 the force dependence of the drift velocity for different temperatures is shown 
m Fig. [TJ In the range of forces of our investigation one can observe two regimes: a low 
force regime (I) in which the velocity increases linear with increasing force and a high force 
regime (II) which exhibit a nonlinear growth of velocity. 

In the first regime, the velocity fulfills the linear response relation 

v = D f3F (2) 

(dashed lines in Fig. [T]) in which Dq stands for the one- dimensional diffusion constant in 
equilibrium and (3 = for the inverse thermal energy. It is of particular interest, that 
the size of the linear regime is temperature dependent so that it decreases with decreasing 
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FIG. 1. Drift velocity v of BLJM65 . The dashed lines correspond to the expected linear response 
behavior (see Eq. [2]). The solid line indicates the critical forces F c>v at which the linear response 
regime (I) pass on to the nonlinear regime (II) (see Eq. [4]). 

temperature. Similar to the work of Williams et alP, we quantify this observation by fitting 
the curves with a symmetric power law 



a 2 F + a . 



(3) 



Using this fit one can define a threshold force 



0.1 ao 



a 2 



(4) 



by the criterion that the dynamical response differs more than 10% from the expected 
linear response behavior. The values of these threshold forces are indicated as solid lines in 
Fig. [Tj The critical forces are further discussed in Ref III D 



To analyze possible finite-size effects for the velocity we have compared the drift velocity of 
the BLJM65 and BL JM1560 at a temperature of T = 0.475. In the linear regime one observes 
that the velocity of the larger system slightly differs from the small system. According to 
linear response theory this just reflects analogous effects for the equilibrium diffusivity as 
reported in RefJ^. Additional differences of the equilibrium diffusivities can result from 
hydrodynamic effects due to the different geometries. Most interestingly, after superimposing 
the velocity of the BLJM65 systems in the linear regime also the results for the nonlinear 
regimes are rather close. Especially the onset of nonlinear effects seems to be equal for both 
system sizes. Thus, understanding the nonlinearity for BLJM65 is sufficient to unravel the 
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FIG. 2. Drift velocity v for BLJM65 and BLJM1560 at a temperature T = 0.475. v was normalized 
via the particular equilibrium diffusion constants -Do- 



underlying physics also of large systems. In particular, this allows one to use the CTRW 
framework in configuration space for the elucidation of microrheological effects. 



B. Drift velocity in CTRW terms 

In principle, any kind of discretization of the trajectories allows the decomposition of 
dynamical quantities in a spatial and a temporal part. In case of the drift velocity one can 
generally write 

< A *II> ,~ 

This relates the long time drift to the average displacement of the tracer particle along 
the force direction during a single elementary step, (Axy) and the average waiting time (r). 
Because the waiting time cannot be affected by the direction of the applied force the value 
of (r) can only depend on even powers of F. Therefore, the linear response regime of small 
forces has to be related to the force dependence of the spatial part. However, in case of the 
nonlinear regime, it is a priori not clear, whether the spatial or temporal effects dominate 
the dynamical responses. 

The force dependence (Ascy ) is shown in Fig. [3j One observes a linear scaling at low and 
intermediate forces. This behavior can be quantitatively understood by taking into account 
that the equilibrium diffusion constant D can be expressed in CTRW terms^ as 
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FIG. 3. Average displacement of the tracer particle during one MB transition (Accy) as a function 
of the applied force F. The solid line indicate the linear response prediction by Eq. [8] while the 
dashed line includes the theoretical prediction of nonlinear effects (Eq. Wj) 



The length scale ao refers to the one-dimensional averaged diffusive length a single A 
type particle moves during one MB transition in an equilibrium system. It is denned by the 
single particle displacement x in one dimension after a large number of transitions n via 

= Um <(*(„) -*(oOT (7) 

ra— s>oo n 

As it was discussed in Ref P^, <2q is a temperature independent quantity. Inserting Eq. [H] 
and Eq. [6] in Eq. [2] one obtains for the linear response regime 

(A*||> = (8) 
which is indicated solid line in Fig. [3] 

For the evolution of (Ash) in the nonlinear regime one can make a prediction by consid- 
ering a one-dimensional periodic potential with a constant distance ao between two adjacent 
minima. Under application of the force, the jumping directions become biased so that one 
gets for the average displacement 



(Axy) = a tanh —j3F. 
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(9) 



10° 



10 1 10 2 

I 



10 3 



FIG. 4. Distribution of MB waiting times r during one MB transition at a temperature T = 0.475. 

This particular ansatz was further discussed by Jack et alP^. By linear expansion, Eq. [9] 
also includes the linear response relation (Eq. [8]). However, from Eq. [9] one would expect 
a negative nonlinear response because the F 3 term of the expansion has a negative sign. 
In contrast, our numerical results display a slightly positive course. One can qualitatively 
understand this behavior by taking into account that in the theoretical ansatz ao is regarded 
as a force independent quantity. Especially at higher forces, one could expect that the 
diffusive length increases with increasing force as well (see also below). Unfortunately, 
one cannot immediately identify a specific length which displays the expected behavior to 
quantitatively reproduce the nonlinear growth of the particle displacement. Thus, we see 
that the spatial aspects of the hopping behavior at large forces cannot be described by this 
simple approach of equidistant minima. 

In any event, the key conclusion from this analysis is the smallness of the nonlinear 
increase. Even at the lowest temperature and highest force it is smaller than a factor of 
1.25. As a conclusion, the spatial part of the drift velocity comes with a distinct linear and 
a very weak nonlinear response. 

Concerning the temporal part we show the waiting time distribution (p(r) i n Fig. |4j In 
the case of small forces, the distribution of waiting times does not change so that the curves 
collapse as expected for the linear response regime. Indeed, in case of higher forces, long 
waiting times do not occur any more so that the distribution functions decay faster. This has 
major consequences for the average waiting time (r) (see Fig. |5|: While for small forces (r) 
stays constant, one can observe a drop for forces higher than a critical force. We characterize 
the critical behavior of the temporal part by determining the critical force of the inverse 
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FIG. 5. Average waiting times (r) as a function of external forces F. 



average waiting time 1/ (r) which is the relevant contribution to the drift velocity (see Eq. [5]). 
It turns out, that the critical force F j_ is identical to F cv . This fact underlines the major 
influence of the temporal part to the nonlinear behavior of the drift velocity. For further 



discussion in |III D| we will thus not distinguish between these critical forces and restrict our 
analysis on F i . 



The transition to the nonlinear regime of v can be regarded as a direct consequence of 
the nonlinear behavior of (r). The observations of the different force dependencies of the 
spatial and the temporal part interestingly illustrate that the discretization of the trajectory 
in the MB approach also leads to a decomposition of linear and nonlinear contributions to 
the drift velocity: In case of small forces, (r) stays constant so that the linear response can 
be completely related to the spatial part. The nonlinear regime, however, is governed by a 
decay of the average waiting time. 



On a qualitative level this picture agrees with the idea of linear response theory as applied 
to a simple ID cosine-potential. Naturally, (Axii) is proportional to the force in the linear 
regime. Thus, any variation of the waiting times, e.g. to a modification of the barriers, gives 
rise to higher-order effects in the resulting mobility. Stated differently, the nonlinearity 
mainly reflects the modification of the underlying PEL upon application of a large external 
force. 
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FIG. 6. Diffusive length scales a?, (two upper lines) and a\ (two lower lines) of the tracer particle 
parallel and perpendicular to the force direction as a function of the external force F. 

C. Diffusion in CTRW terms 



As it was already mentioned above (see Eq. |6J), the one-dimensional diffusion constant 
in equilibrium of a single particle is given in CTRW terms as the ratio of the isotropic 
diffusive length a 2 , and the average waiting time (r). Going towards driven system one has 
to distinguish between the diffusive processes parallel and perpendicular to the force direction 
which were characterized by the diffusion constants D\\ and D±, respectively. Because (r) is 
a universal quantity, expected differences between parallel and perpendicular diffusion are 
related to the anisotropy of the respective length scales. In analogy to Eq. [7| we define the 
length scale in perpendicular direction a\ as 



lim 



((x ± (n)-x ± (0)) 2 ) 



n 



(10) 



and the parallel length scale a 2 as 



lim 

n— >oo 



((x||(n)-xn(0)) 2 )-(( 3 ;||H-X||(0))r 
n 



(11) 



In the latter definition we subtracted the systematic shift of the tracer particle along the 
force direction. The corresponding length scales are shown in Fig. [6] One can observe that 
both, a\ and a», exhibit an increase with increasing force which is stronger for the parallel 
direction. Interestingly, although a\ and a| are temperature independent in equilibrium, 
corresponding to the linear regime, the degree of nonlinearity strongly depends on temper- 
ature. As already indicated above one may expect that the PEL properties change in the 
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FIG. 7. Nonlinear responses of the diffusion coefficients parallel (D\\) and perpendicular (D±) to 
the force direction and of the drift velocity v at a temperature T = 0.475. 

nonlinear regime. Since these modifications are expected to be strongly anisotropic it is not 
surprising that observables, defined parallel or orthogonal to the applied force, may behave 
differently. 

Dividing a]_ and ajj by 2(r) one obtains the parallel and the orthogonal diffusion constants 
Dm and D±, respectively. Since the denominator is always the same the different increase 
of a| as compared to a\ directly translates into the corresponding relative increase of Dy 
as compared to D±. Furthermore, the weak force-dependence of (Axy) directly translates 
into a weaker force dependence of the mobility as compared to both diffusion constants. 
Actually, the observed relation Dy 3> D± > jp (see Fig. [7j) was already reported for the 
lithium dynamics in a lithium silicate systempH. This suggests that the force dependent 
anisotropy of the relevant length scales is a general property of disordered systems. 



D. Critical forces 

Finally, it is of interest to compare the nonlinear response of the different observables on 
a more quantitative level. Since we are mainly interested in the onset of nonlinearity an 
appropriate measure is the critical force as introduced in Eq. 4 Here we concentrate on o|, 
a\ and 1/ (r) as the key observables of the CTRW approach. The remaining variable (Ax\\) 
is omitted because its nonlinearity is extremely small. Furthermore, the critical force of the 
drift velocity F CfV is not discussed separately because, as it was mentioned above, it has been 
found to be identical with F i . 
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FIG. 8. Critical forces F C) x of different dynamical quantities X = ajj, and ^. The dashed lines 
indicate a behavior F c x oc T L5 or rather F c x oc T 2 5 . Inset: F i as a function of ;y S^ with the 
equilibrium diffusion constant Dq. The dashed line indicates a linear behavior. 

In Fig. [8] all critical forces are shown as a function of temperature. It can be seen that in 
general for a given temperature each observable has a different critical force. Furthermore 
the temperature dependence varies quite significantly. Whereas for the length scales the 
critical forces scale similarly to F c> x oc T 15 , the critical force of the inverse waiting time 
roughly scales like F c i oc T 25 . 

A theoretical prediction about the onset of nonlinearity can be found in^ 2 . It was sug- 
gested by the authors that linear response relation of the occurrent flux in a driven system 
holds exactly when the steady state fluctuation theorem 
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p(v(t) = -A) 
is fulfilled. From this condition the the scaling 



p(v(t)-A) ^ exp{Am (12) 



^c. « ^ (13) 



has been derived in which Dq stands for the equilibrium diffusion coefficient. In the inset 
of Fig. |8j one can observe that F c i , and therefore naturally F C)V as well, fulfills the scaling 
prediction. 
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IV. CONCLUSION 



In the present paper we have analyzed the dynamics of a driven single particle in a 
supercooled liquid in terms of a CTRW. In case of the drift velocity the discretization of 
the MB approach allows one to identify the temporal part of the CTRW as the key source 
of nonlinear response, whereas the spatial part only exhibits a weak positive contribution. 
This positive contribution can be qualitatively understood by the fact that the distances 
between adjacent MB increase due to the application of higher forces. Since the finite-size 
effects are small these results also reflect the origin of the nonlinear response in the much 
larger system. 

In terms of a CTRW analysis we have been able to determine the relevant diffusive lengths 
of the tracer particle diffusion along the directions parallel and perpendicular to the force. 
Among others, this enables the definition of a parallel diffusion constant which s not directly 
accessible from real space trajectories^. 

The force dependencies of parallel diffusion, perpendicular diffusion and drift velocity are 
significantly different. Thus, there is no general critical force which designates a mutual 
transition to a nonlinear regime. However, in terms of dynamical quantities, one can regard 
the critical force of the inverse waiting time as a key entity because of its relation to the 
velocity (neglecting the minor nonlinear effects of the spatial part). For this quantity one 
also observes the validity of the theoretical scaling prediction by Evans et al.. 

This observation is of major conceptual interest because the theoretical prediction is 
based on steady state fluctuation arguments while the CTRW waiting time is determined by 
the underlying PEL, i.e. by the distribution of minimum energies and saddle heights as well 
as the topography. Because the waiting time decreases under the application of the external 
force, one can assume that the force dependence of at least one of these distributions is 
responsible for this behavior. It remains to be shown which specific variations of the PEL 
give to the transition for linear to nonlinear dynamics. 
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